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AC Josephson effect in the long voltage-biased SINIS junction 
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Theory of non-stationary coherent effects is developed for superconductor-normal-superconductor (SNS) 
structures with relatively strong normal scattering on S/N interfaces (interface resistance is large compared 
to intrinsic resistance of N metal). Analytical expressions are found for the time-dependent anomalous Green 
functions induced in the N region under the fixed-voltage-bias. The amplitude of the current oscillations is 
determined in non-equilibrium conditions. Non-stationary correction to the distribution function is calculated 
in high-temperature limit and found to be slowly decreasing with the temperature, leading to the dominance of 
the second-harmonic term in the Josephson current, I s {t) oc sin(4eVt) at high temperatures and low voltage. 
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1. Coherence effects and general equations. 

The superconducting hybrid (normal metal - supercon- 
ductor) structures are very rich systems that have been 
studied for the past few decades both theoretically and 
experimentally (cf. [T] for a relatively recent brief re- 
view). The proximity effect has been shown to induce 
superconductive correlations into normal part of the sys- 
tem, where they decay over the length £ e = ^hD/e 
that can be large at low energies. One of the conse- 
quence of the proximity effect is the ability of the nor- 
mal metal to carry out phase-sensitive current if length 
L of the junction is smaller then £ eo at characteris- 
tic energies erj. Equilibrium Josephson effect in diffu- 
'sive SNS systems is well understood by now and can 
be well explained in terms of stationary Andreev levels 
or more directly through Green's functions approach. 
When the constant voltage is applied to the junction, 
nonstationary Josephson effect arises. If the voltage is 
high, V ^> Etii, than effects of coherence between su- 
perconducting reservoirs can be neglected and current is 
stationary with the rich subgap structure [5] • When the 
voltage is not very high, time-dependent contribution to 
the current can become important, even at rather high 
temperatures [3j |H [6] . 

In this paper we consider nonstationary Josephson 
effect in a long symmetric voltage-biased SNS junction 
with large interface resistance, so that r = Rb/Rn ^ 1 
(here Rb is the single barrier resistance and Rn the 
resistance of the normal wire). This problem was first 
considered in [7] , where Josephson current up to the first 
order in r _1 was calculated. It was shown, that Joseph- 
son current decays ex ponentially when temperature is 
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high comparative to the inverse diffusion time. How- 
ever, it was shown later [3j, that due to nonequilibrium 
effects time-dependent current does decay as slow as 
T _1 at high temperatures. Here we microscopically cal- 
culate the current, with all nonequilibrium effects taken 
into account. In doing so, we have to include the terms 
which are formally of higher orders in r ; however the 
resulting current is not necessary smaller than lowest- 
order term, but can dominate it, as explained below. 

To describe this system we use Keldysh method, de- 
veloped for superconductivity by Larkin and Ovchin- 
nikov 8 . In this method the Green function is 4 x 4 
matrix in Keldysh and particle/hole 2x2 spaces. This 
matrix G contains all the information about spectrum of 
the system and on the distribution of electrons over the 
energy levels. The resulting Usadel equation is written 
in terms of the disorder-averaged semiclassical Green 
function G = G(ii,i2> r ) (f° r the detailed review, see 
[9]). This equation is presented below, it can be used 
to describe any nonstationary phenomena with low en- 
ergy scales involved (compared to Fermi energy). It 
involves time convolution operation: (/ o g) t 2 ) = 
J_ f(ti,t)g(t,t2)dt. It's convenient to introduce t = 
(t\ + t2)/2, t = t\ — t%. Then Usadel equation in the 
absence of the vector potential (which is supposed to be 
zero throughout the paper) and electron pairing reads 
as follows (we work in units e = &B = /i=l): 



Dd x {God x G)+d T [a 3 ,G] 



G} + icp. 



G = I 
(1) 

with er 3 = If 3 and I = —i (£j n o G — G o Sj ra J where 
Si„ is self-energy matrix, that accounts for inelastic 
and dephasing processes. Time-dependent electric po- 
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tential (p(t) enters this equation through (ti, ^2) = 
tpfe) — <p(ti). It has to be determined self-consistently 
via electroneutrality condition tp(t) = jtrG (t,t) [9]. 

Keldysh Green function G obeys normalization con- 
dition G o G — lS(ti — t^} which allows for the anzatz 
G K = G R o h — h o G A where h stays for the matrix 
distribution function. The latter can be chosen to be di- 
agonal: h = hof + /13T 3 . It is convenient to implement 
Fourier transform over r, i.e. to pass to Wigner repre- 
sentation. This way, we write G(t, t) = J G(e,t)e~ leT de. 
We note that the convolution operator simplifies in the 
mixed representation due to the following identity: for 



/(e), g(e,t) = e-^ t g(e), one gets 



(fog) (M) = 



y)ff(e- y) 



Our goal is to find G R , G A ,G K for the SINIS struc- 
ture and to calculate the current. Current density is ex- 
pressed through Keldysh component of the matrix cur- 
rent j — G o VG as follows: 



■tTT 6 j K {t,t) 



(2) 



Matrix equation ([I} has to be supplemented with the 
Kupriyanov-Lukichev [10] boundary condition at both 
SIN interfaces. At the right interface it takes the form: 



2RsN<JNj = [Go, G r 



(3) 



with Rsn being the interface resistance of the barrier 
per unit area in the normal state. We use indices I, r for 
the left and right reservoirs. 

Now we have to work out retarded, advanced and 
Keldysh components of equations (| 1 13j) . In the explicit 
form, they are presented in the Appendix. Retarded and 
advanced components determine the generalized spec- 
tral properties of the system and Keldysh components 
describe distribution of electrons over this generalized 
spectrum. In what follows, we simplify these equations 
using the smallness of the parameter r -1 , as described 
in the next section. 

We consider the constant-voltage-biased setup and 
neglect the spatial dependencies of all the relevant 
quantities in the directions, perpendicular to the wire. 
We measure the length in units of L, assuming that 
bulk superconductors are situated at points x = — \ 
and 1 = i at the voltages — ^ & n d -j correspond- 
ingly. We suppose that bulk superconductors are good 
reservoirs at temperature Ts] under this assumption 
their Green functions Gi_ r can be obtained from 
the standard BCS form Gbcs by the gauge trans- 
formation Gi, r (tl,t 2 ) = SLr{tl)GBCs(tl,t2)Sf r (t 2 ) 



with Gg^, — 



= (g s r a + fsr 1 



R(A) 



U BCS 



T BCS 



and Si <r (t) 



tanh(_y (Gi cs - G 
traigh 

g R (e T V/2) e* ivt f R (e 



exp(± 1 |^f 3 )i. Straightforward calculation gives: 



<3f r (f,e)= I "* v : ' " I (4) 

^ ' \ e±* vt f R (e) -g*(e±V/2) » 



and 



tanh( 



e±V/2s , 
2T S > ' 

tanh(- 



±V/2 \ 
2T S > 



(5) 



The explicit form of the energy-dependent functions 

gf A \ ff A) is: 



gf A) (e) = - (± VS - z£ s ) , fs^'(e) = ± ivs, (6) 



MA), 



where 

A sign e . ., . . . . 
Vs = 7 =JL=6(\e\-A), 



VA 2 - e 2 



9(A-\e\) 
(7) 

2. Spectral functions. Here we calculate Green func- 
tions G^ R,A ^ which describe proximity effect in normal 
region. We write matrix Green function G RtyA ^ (x, e, t) 
in the form 



qR(A) 




R{A)s 



R(A) 



f 

f? W Td-g? {A) ) 



(8) 



Since the proximity effect is weak due to the presence 
of barriers, we can expand the equations for G^ R,A ^ in 
powers of small parameter r _1 . To the first order one 
finds: 

9?,2~r~ 2 , f« 2 A ~r-i (9) 

Besides, we assume here that (p— ~ r~ 2 , which is self- 
consistent assumption. Smallness of <p is due to the fact 
that almost all voltage drops at the barriers at r 1. 
We consider Eqs. (|35l36p and keep terms, proportional 
to r _1 in the equations, and terms of the order of unity 
in the boundary conditions for anomalous Green func- 
tion f R (x,e,t), to obtain: 



±r -l e ±iVtfR 



Here 



2l€ 



Th 



(10) 



(11) 



and 7 = (TinExh) -1 is the dimensionless inelastic scat- 
tering rate. Linear approximation leading to Eqs. (|10p is 
applicable at all energies e, if (-fr) -1 ~ r in E g <C 1 (here 
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Eg is the proximity-induced minigap [11] at 7 = 0). 
Other Green functions are expressed via relations 



ff a {x,e,t) = fMx,-€,t) 



(12) 



f^ A \x,e,t) = f^>(x,e,-t) 



R(A), 



The above equations are easy to solve, with the follow- 
ing result: 



fx 



iVt 

e cos k £ 



, l\ , -iVt f 1 

x + - J + e cos k e I x - - 



(13) 



For the future convenience, we introduce definitions 



u 



_£iM u ( e ) i V R 



(14) 



^th u(e) = v(e) = , so that 

fi^\ x = \)= u R We ±iVt + v R We* iVt . (15) 

Retarded and advanced normal Green functions can be 
found from the normalization condition: 

<?f = \f? off, 9? = \f?of« (16) 
Electric potential is 

tp(t) = ~ J {h 3 + g R oh -h o g$) de (17) 

with g R = i [f R o, /^j . Linearization over f\ 2 breaks 
down at energies too close to the gap A: \e — A| < 
^r- max(£ , 7'? l /A, 1), but the resulting square-root sin- 
gularity |e — A| -1 / 2 does not influence our results. 

3. Kinetic equation and its solution. In this 
section we determine the electron distribution function. 
Its form is determined by interplay of several processes. 
The electrons diffuse from one superconducting lead 
to another during characteristic time td = and 
are subjected to normal and Andreev scattering at the 
boundaries. Electron-electron scattering tends to ther- 
malize them to equilibrium with some effective temper- 
ature T e . Electron-phonon scattering tends to bring T e 
closer to the substrate temperature Ts, thus taking the 
energy out of electron system. Another (and more effec- 
tive at low temperatures) channel of electron cooling is 
the tunneling of hot electrons to superconducting reser- 
voirs [12]. In the SINIS junction, the role of inelastic 
scattering is relatively large, since Andreev reflections 
are suppressed due to weakness of the proximity effect, 
so that electron spends a long time, bouncing back and 



forth between the boundaries. At % = jr 2 ^S> 1 the 
electron distribution thcrmalizes and reads 



h(x,e,t)^h { °\e)^ + O( X - 1 ) 



(18) 



where h^\e) — tanh(7^) is equilibrium distribution 
with some effective temperature T e ; in general, T e ^ T$. 
Effective temperature T e has to be determined from the 
heat balance equation (cf. e.g. [13]). After that is done, 
one can calculate non-equilibrium correction to the dis- 
tribution function (the second term in (fl8)) '). Note, that 
it can be important in terms of calculating the Joseph- 
son current, even if it small. The point is that thermal 
distribution leads (see below) to the amplitude |/ s | of 
the ac Josephson current I(t), which is exponentially 
small in £/£t- O n the other hand, non-equilibrium 
corrections decay with temperature and length much 
slower. These non-equilibrium corrections are rather 
interesting, since they arise due to coherent Andreev 
reflections. 

In order to find this nonequilibrium correction, we 
simplify general kinetic equations (|37[) . (|39[) by separat- 
ing terms of different orders over small parameter r , 
and adopting the simplified form of collision integral, i.e. 
" T-approximation" . In doing so, one has to consider the 
boundary condition (|39|) up to the r _1 terms and kinetic 
equation (|3"T|) up to the r~ 2 terms. Besides, these are 
small energies e <C A, where deviations from the equi- 
librium are important, so we put 775 (e) = 0, <ts( e ) = 1- 
We look for the distribution function in the following 
form: 



h (x,e,t)=h ( ° ) (e)+r- 2 h ( o\e,x,t) 
h 3 (x,e,t) = r- 2 h { ^\e,x,t) 



(19) 



(2) (2) 

Next, due to the spatial symmetry, h 3 and are 
odd and even functions of x, correspondingly. Taking 
this into account, we write boundary conditions (|39|) 
only at the x = \ boundary, to obtain with the neces- 
sary accuracy: 



40*fc$l*=i = J i T J2 



(20) 



with J12 terms which depend on only, and are 
known therefore: 

Jifi = r (A fl 2 o e* vt X + e ±lVt X o /£) (21) 

with X(e) = i (4+ - 

We used here special notation for the " energy-shift" 
subscript: for any function of energy /(e) we define 

/±(e) = /(e±^) f++ = f(e + V) f__=f( e -V). 

(22) 
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Since the boundary conditions are known, we proceed 
with the kinetic equation. Weakness of the proximity 
effect allows to neglect modifications of the diffusion co- 
efficient, and also the terms, mixing ho and /13 in the 
matrix equation (|37|) . Adopting r-approximation for 
the collision integral, we obtain 



E Th d x h 



21,(2) 



(d4 2) +^4 2) ) = o 



r _1 at low temperatures and exponentially suppressed 
at high temperatures. In the limit of large supercon- 
ducting gap A 3> max(-Br/i) V) the result reads 



de 



with 



lt(t) = (Ki,2°e* iVt 



e TlVt o K h2 ) 



(26) 



(27) 



E Th dlh^ - (d t h^ + + ir^-h^) = 

(23) 

Kinetic equations (j23|) are to be written for each har- 
monic (in terms of the total time t) of the distribution 

function, with the electric potential term that couples 
(2) (2) 

h and h 3 due to the self-consistency condition (see 
(fT7| ). In the limit of large temperatures T e ^> Eth, 
when the nonequilibrium contribution to the current 
becomes important, this term is exponentially small in 
L/£t and we neglect it. Besides, if V <C T e , one has 
X = V/AT e . 

We look for the nonequilibrium correction in the fol- 
lowing form: 



h% (x, e, t) = J2 n =Q ±1 cos(K nV x)e 



hf\x, e, t) = 2 n=0 ±1 #n, e sin(K„yx) 



-2iVnt 



-2iVnt 



(24) 



Boundary conditions (|2TJ)) allow to calculate A, B: 



A - l v A 



B,, 



16 T e K^y sin^Knv) 
r V B n , e 

16 T e K n v COs(^K Tl v f ) 



(25) 



with 



A(n, e) = (a+ - a_ ) 5„, + (w+ - u*) - (u_ - «+) <5 n ,-i 

B{n, e) = (a+ + a_) 5 n> o + («+ + «-) S n ,i + («- + «+) <5n,-i 

where a = u fl + u A . 

4. Time-dependent current 

With the distribution functions (fT5|) being deter- 
mined, we turn to the calculation of time-dependent 
electrical current. It is convenient to calculate it in the 
vicinity of the right boundary, with j K in There 
are lot of terms in @ , but we keep only those of them, 
that are proportional to the first order contribution to 
the anomalous function G RtyA ^ and contribute to the os- 
cillating part of the current. This allows us to get the 
leading terms at both low and high effective tempera- 
tures comparative to Exh- More precise treatment will 
give some corrections, that are of the higher orders of 



with R being the resistance of the SINIS junction in 
the normal state: B = 2Bb- We have introduced here 
K\2 = fi2 o h 2 — hi o f A 2 : where anomalous and distri- 
bution functions are supposed to be calculated at x = h. 
To proceed, we note that I(t) it can be considered as a 
sum of two different contributions, according to fl8|. 



Equilibrium current. It reads: I eq = 5R[e 



-2iVt 



h] 



with 

h = 



1 

2Br 



[v(e + V/2) - vi-e + V/2)] tanh —de 

6 (28) 

We note that at zero-voltage limit one has 
= tanh(2^-) and for the amplitude I\ we get 
usual result 7 for the equilibrium critical cur- 
rent: Ix(V — 0) = —gz. f hs^sv de. In general, 
the integral (|28p can be reduced to the sum over 

2 " T = 1 where 



residues to give 1\ = — - 



Rr 



n=l 



q n sinh q n 



2(2n-l)-!rT e -iV . 

g" = V E Th +7- 



In the limit of T e 3> Eth and 7 < 1 we obtain 



\h\=c- 



b^ 



2-kT, 



2 Ei 



4tt, V < T e 



At very low temperatures tanh ^4r 
integral for the current can be done explicitly 



(29) 

signe and the 



\h\ 



IE' 



Th 



rB 



In cot 




(30) 



The above expression is valid for any voltage if Trj„ <C 1 
and for high voltages Wi n ^> 1 otherwise. 

Nonequilibrium current. Calculation starting 
from Eq. (|2"6]) leads to: 



L neq 



-2iVt 



-4iVt; 



with complex amplitudes $1,2 
4>\ 2 are given by: 



£ <f 2 ] (31) 
■^^01,2, where 



0i = r (xo - Do) - rv(io + y ) + (xv + yv)(T 2 v - Tv) 



= - (xv + yv)T2v 



(32) 
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with x t 



- cot( ^; /2) , Ve = tan( K K ; /2) and r e = 
; 1 — ; The results (l32t are ap- 

s/2(ie/E Th --y)sin- s /2(ie/E Th - 1 ) 

plicable for any values of V/Exh ratio, but in the lowest 
order over V/T e . 

Here we analyze the case of weak inelastic scattering 
7< 1. In the limit of zero voltage, |^| = 27 _1 |0i| ^> 
|0i | and the second harmonic dominates. For small volt- 
ages V <SC E Th one has: 



I'M 



$l = - 



and 



irE Th V 
WRr 3 T e 



nE Th V 
16i?r 3 T e 



+ 4j- 4 {V/E Th ) 2 



< V 



(33) 



\{E Th /Vf 



n n « v 

(34) 

Comparing first lines of Eqs. ([34f and ([29]) we find that 
non-equilibrium second harmonics dominates the ac 
current at T e > Ej-h hi 2 (fr) and low voltages V <C r^ 1 ■ 
Similar phenomenon was observed in Refs. |14[ 115] . 
where subharmonic Shapiro steps with slowly decreasing 
(upon T increasing) amplitudes were found. Qualita- 
tive theory of this phenomenon was proposed by Arga- 
man [3], who discussed it in terms of time-dependent 
Andrccv bound states with non-equilibrium popula- 
tion. Our result (|34|) contains the same V/T e depen- 
dence at low V as found in Ref. g]; however, we got 
Ineq ~ X -1 -JiJef aman '• We expect therefore that the re- 
sult of [4] is valid (up to numerical factor of order unity) 
under the condition \ < 1, which we do not consider 
here. To understand the origin of the whole effect it 
is useful to note that the result [4] for the second har- 
monics is very similar to (9) of Ref. [6] where Debye 
relaxation contribution to the dc conductance of SNS 
junction was estimated. Moreover, preprint version of 
Ref. [6] contains the same kind of estimation for the 
SINIS junction, with the result oc rf n , very much like 
our expression for $2- We believe therefore that non- 
equilibrium second harmonics of the current originates 
from the same Debye relaxation mechanism. 

Non-equilibrium ac current beyond linear in V ap- 
proximation never was calculated previously, to the best 
of our knowledge. Eqs. (|33l34p demonstrate that at not 
very low voltages non-equilibrium first harmonics of the 
current <&i becomes comparable to $2 and then exceeds 
it. Note also 7r-shift in phases of $12 with growth 
of voltage. Technically, $1 originates from the peaks 

in nonequilibrium stationary parts of distribution func- 

(2) 

tions, /193(e) at e = ±V/2. These peaks result from 



0.5 y- 



Nonequilibrium amplitudes |0i(V)| (normal line) and 
|02(V)| (dashed line) are calculated with Eq. (|32[) for 
7 = 0.2 and V < T e . 



the modulation of the spectral density with the Joseph- 
son frequency, cf. Eq. (fl~5|) . For the particular value of 
7 = 0.2, the absolute values of amplitudes 0i,a(V) at 
low voltages V < Etk "C T e are plotted in the Figure. 
5. Conclusions 

In this paper we have considered the AC Joseph- 
son effect in a long SINIS junction in the case of tem- 
perature and voltage low with respect to the bulk gap 
A. The main results are given in Eqs. (|32l33l34j) for 
non-equilibrium contribution to the current in the high- 
temperature range T e 3> Eth- At high electron temper- 
atures T e in the normal wire, the nonequilibrium con- 
tribution to the current becomes dominant. We found, 
that the ac current contains two harmonics: the first 
with basic Josephson frequency uij — 2eV/H and the 
second harmonic 2lo j. Second harmonic has the largest 
amplitude at low voltage and high temperature. Note, 
that calculation of non-equilibrium effects at low tem- 
peratures is complicated due to the necessity to account 
for the time-dependent electric potential in the normal 
wire; we leave this problem for the future studies. 

We are grateful to Ya. V. Fominov for many use- 
ful advises, and to H. Bouchiat, H. Courtois, S. Gueron 
and V. V. Ryazanov for illuminating discussions. This 
research was supported by RFBR grant 07-02-00310 and 
by the RAS Program "Quantum physics of condensed 
matter" . 

Appendix. 

Here we explicitly write out nonstationary Usadel 
equation and boundary conditions. Retarded compo- 
nent of (JTJ) takes the following form: 



^Dd x (G R od x G R )+d T [r 3 ,G R } + ^d T {r\G R }+i^G R = I R 



(35) 
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with I R = —i 



£* o, G R 



It is supplemented with re- 



tarded component of boundary condition ([3]) : 

2R SN a N tf = [Gf-o,G R ]. (36) 

Since advanced components are identical (with R — > A 
substitution), we proceed with kinetic equation, which 
results from the Keldysh component of Usadel equa- 
tion: 

—D \d x (d x h -G R ° dji o G A ) + (j R od x h~d x ho j A )] 

+ (df o drh + iG R o (p_h) - (d T h oGf+ iip-h o G A ^j = I st 

(37) 

with 



G R(A) = \ {G^),f 3 } , I st = -i(g r ocj-&o G a 



& = Y> R oh-hoY 1 A -Y, K . 

(38) 

The Keldysh component of boundary conditions ^ 
reads 



"2Rsn<J N (djii - Gf o d x hi o Gf ) = (G R ou~uoG a 

(39) 

where u = G R o Sh — Sh o G A , Sh = h r — hi 
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